pest injuires network

To record the the co-occurence network of pest injuries by country

library(ggplot2)
require(dplyr)

## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## 
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## 
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

require(reshape)

## Loading required package: reshape
## 
## Attaching package: 'reshape'
## 
## The following object is masked from 'package:dplyr':
## 
##     rename

require(reshape2)

## Loading required package: reshape2
## 
## Attaching package: 'reshape2'
## 
## The following objects are masked from 'package:reshape':
## 
##     colsplit, melt, recast

require(qgraph)

## Loading required package: qgraph

library(gridExtra)
library(cowplot)

## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:ggplot2':
## 
##     ggsave

library(lubridate)

## 
## Attaching package: 'lubridate'
## 
## The following object is masked from 'package:reshape':
## 
##     stamp

library(XLConnect)

## Loading required package: XLConnectJars
## XLConnect 0.2-11 by Mirai Solutions GmbH [aut],
##   Martin Studer [cre],
##   The Apache Software Foundation [ctb, cph] (Apache POI, Apache Commons
##     Codec),
##   Stephen Colebourne [ctb, cph] (Joda-Time Java library)
## http://www.mirai-solutions.com ,
## http://miraisolutions.wordpress.com

You can also embed plots, for example:

Filepath <- "E:/Google Drive/1.SKEP1/SKEP1survey.xls" # please check your file in shared google drive
Filepath <- "~/Google Drive/1.SKEP1/SKEP1survey.xls"
data <- readWorksheetFromFile(Filepath, sheet = 1)

#### clean define the missing value ####
data[data == "-"] <- NA # replace '-' with NA
data[data == ""] <- NA # replace 'missing data' with NA
#### end cleaning of data ####

#### to lower variable names ####
names(data) <- tolower(names(data))
#### end setting the varibales ####

data <- transform(data, 
                  phase = as.factor(phase),
                  fno = as.character(fno),
                  identifier = as.character(identifier),
                  country = as.factor(country),
                  year = as.factor(year),
                  season  = as.factor(season),   
                  lat = as.numeric(lat),
                  long = as.numeric(long),      
                  village = as.character(village), 
                  fa = as.numeric(fa),
                  fn = as.character(fn),
                  lfm = as.character(lfm),
                  pc = as.factor(pc),
                  fp = as.character(fp),        
                  cem = as.factor(cem),     
                  ast = as.factor(ast),       
                  nplsqm = as.numeric(nplsqm),
                  ced = dmy(ced),# Date data try to use as.Data(., format = '%d-%b-%y') it is not working
                  cedjul = as.numeric(cedjul),
                  hd = dmy(hd), 
                  hdjul = as.numeric(hdjul),     
                  ccd = as.numeric(ccd),
                  cvr = as.character(cvr),
                  vartype = as.factor(vartype),
                  varcoded = as.factor(varcoded),
                  fym = as.character(fym),
                  fymcoded = as.factor(fymcoded),
                  n = as.numeric(n),
                  p = as.numeric(p) ,
                  k = as.numeric(k),
                  mf = as.numeric(mf),        
                  wcp = as.factor(wcp),      
                  mu = as.character(mu) ,     
                  iu = as.numeric(iu),     
                  hu = as.numeric(hu),      
                  fu = as.numeric(fu),      
                  cs  = as.factor(cs),      
                  ldg  =  as.numeric(ldg),  
                  yield = as.numeric(yield) ,
                  dscum = as.factor(dscum),   
                  wecum = as.factor(wecum),   
                  ntmax = as.numeric(ntmax), 
                  npmax = as.numeric(npmax),    
                  nltmax = as.numeric(nltmax),  
                  nlhmax = as.numeric(nltmax),  
                  waa = as.numeric(waa),      
                  wba = as.numeric(wba) ,   
                  dhx =  as.numeric(dhx),  
                  whx =  as.numeric(whx),     
                  ssx  = as.numeric(ssx),  
                  wma = as.numeric(wma), 
                  lfa = as.numeric(lfa),
                  lma = as.numeric(lma),   
                  rha  = as.numeric(rha) ,
                  thrx = as.numeric(thrx),    
                  pmx = as.numeric(pmx),    
                  defa  = as.numeric(defa) ,
                  bphx = as.numeric(bphx),   
                  wbpx = as.numeric(wbpx),    
                  awx  = as.numeric(awx), 
                  rbx =as.numeric(rbx),   
                  rbbx = as.numeric(rbbx),  
                  glhx  = as.numeric(glhx), 
                  stbx=as.numeric(stbx),    
                  rbpx = as.numeric(rbpx), 
                  hbx= as.numeric(hbx),
                  bbx = as.numeric(bbx),    
                  blba = as.numeric(blba),    
                  lba = as.numeric(lba),    
                  bsa = as.numeric(bsa),    
                  blsa = as.numeric(blsa),  
                  nbsa = as.numeric(nbsa),  
                  rsa  = as.numeric(rsa),   
                  lsa = as.numeric(lsa),    
                  shbx = as.numeric(shbx) ,  
                  shrx = as.numeric(shrx),    
                  srx= as.numeric(srx),    
                  fsmx = as.numeric(fsmx),   
                  nbx =  as.numeric(nbx),   
                  dpx = as.numeric(dpx),    
                  rtdx  = as.numeric(rtdx),  
                  rsdx  = as.numeric(rsdx),
                  gsdx  =as.numeric(gsdx),   
                  rtx = as.numeric(rtx)
) 

## Warning: All formats failed to parse. No formats found.

## Warning: NAs introduced by coercion

#### Delete the unnessary variables variables without data (NA) ####
data$phase <- NULL # there is only one type yype of phase in the survey
data$identifier <- NULL # this variable is not included in the analysis
data$village <- NULL
data$fa <- NULL # field area is not include in the analysis
data$fn <- NULL # farmer name can not be included in this survey analysis
data$fp <- NULL # I do not know what is fp
data$lfm <- NULL # there is only one type of land form in this survey
data$ced <- NULL # Date data can not be included in the network analysis
data$cedjul <- NULL
data$hd <- NULL # Date data can not be included in the network analysis
data$hdjul <- NULL
data$cvr <- NULL
data$varcoded <- NULL # I will recode them 
data$fym.coded <- NULL
data$mu <- NULL # no record
data$nplsqm <- NULL
#### Delete the unnessary variables variables without data (NA) ####

OutputProfile <- data %>% 
        select (
        country, 
        season,
        dhx,  
        whx,     
        ssx,  
        wma, 
        lfa,
        lma,   
        rha, 
        thrx,    
        pmx,    
        defa,
        bphx,   
        wbpx,
        awx, 
        rbx,   
        rbbx,  
        glhx, 
        stbx,
        hbx,  
        bbx,    
        blba,    
        lba,    
        bsa,    
        blsa,  
        nbsa,  
        rsa,   
        lsa,    
        shbx,  
        shrx,    
        srx,    
        fsmx,   
        nbx,   
        dpx,    
        rtdx,  
        rsdx,
        gsdx,   
        rtx 
) 

data <- OutputProfile

#### The Philippines n = 40 ####
php <- data %>% 
        filter( country == "PHL") %>%
        select(-c(country, season))

php <- php[,apply(php, 2, var, na.rm = TRUE) != 0] # exclude the column with variation = 0

php <- php[complete.cases(php),] # exclude row which cantain NA
#### end Philippines subset ####

#### India  N = 105 #### 
ind <- data %>% 
        filter( country == "IND") %>%
        select(-c(country, season))

ind <- ind[,apply(ind, 2, var, na.rm = TRUE) != 0] # exclude the column with variation = 0

ind <- ind[complete.cases(ind),] # exclude row which cantain NA
#### end India subset ####

##### Indonesia N= 100 ####
idn <- data %>% 
        filter( country == "IDN") %>%
        select(-c(country, season))

idn <- idn[,apply(idn, 2, var, na.rm = TRUE) != 0] # exclude the column with variation = 0

idn <- idn[complete.cases(idn),] # exclude row which cantain NA
#### end Indonesia subset ####

#### Thailand n = 105 ####

tha <- data %>% 
        filter( country == "THA") %>%
        select(-c(country, season))

tha <- tha[,apply(tha, 2, var, na.rm = TRUE) != 0] # exclude the column with variation = 0

tha <- tha[complete.cases(tha),] # exclude row which cantain NA
#### end Thailand subset ####

#### Vietnam n = 105 ####
vnm <- data %>% 
        filter( country == "VNM") %>%
        select(-c(country, season))

vnm <- vnm[,apply(vnm, 2, var, na.rm = TRUE) != 0] # exclude the column with variation = 0

vnm <- vnm[complete.cases(vnm),]

abbre <- c("DH", "WH", "SS", "WM", "LF", "DEF", "BPH", "WPH", "AM", "RB", "GLH", "BLB", "LB", "BS", "BLS", "NBS", "RS", "LS", "SHB", "SHR", "SR", "FSM", "NB", "DP", "RTD", "RT")
php$ldg <- NULL
names(php) <-abbre
InsectInjuiry <- c(1:11, 26)
Disease <- 12:25
injuries <- list(InsectInjuiry, Disease)
#==========#
cormat <- cor(php, method = "spearman", use = "pairwise") # spearman correlation
php_qnet <- qgraph(cormat, 
               layout = "spring", 
               threshold = 0.20,
               maximum = 1,
               groups = injuries, 
               color = c("skyblue", "wheat"),
               vsize = 5, 
               line = 3,
               diag = FALSE,
               posCol = "forestgreen",
               negCol = "firebrick3",
               borders = FALSE,
               # legend = TRUE,
               vTrans = 200,
               #nodeNames = Names,
               legend.cex = 0.3,
               title = "Philippines"
               #filetype = "png",
               #filenames = "figs/ind_network.png"
)

qgraph(cormat, 
               layout = "spring", 
               threshold = 0.20,
               maximum = 1,
               groups = injuries, 
               color = c("skyblue", "wheat"),
               vsize = c(1.5,10), 
               line = 3,
               diag = FALSE,
               posCol = "forestgreen",
               negCol = "firebrick3",
               borders = FALSE,
               # legend = TRUE,
               vTrans = 200,
               #nodeNames = Names,
               legend.cex = 0.3,
               title = "Philippines"
               #filetype = "png",
               #filenames = "figs/ind_network.png"
)

# compute the topological properties
cen <- centrality_auto(php_qnet)
topo <- cen$node.centrality
topo$node <- row.names(topo)
row.names(topo) <- NULL
topo.melt <- melt(topo)

## Using node as id variables

p1 <- topo %>% ggplot(aes(x= Strength, y = reorder(node, Strength))) +
         geom_point(size = 3, color ="red") +
        theme_bw() +
        theme(panel.grid.major.x =  element_blank(),
              panel.grid.minor.x =  element_blank(),
              panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
        xlab("Node degree") + 
        ylab("Variables")  

p2 <- topo %>% ggplot(aes(x= Betweenness, y = reorder(node, Betweenness))) + 
        geom_point(size = 3, color ="blue") +
        theme_bw() +
        theme(panel.grid.major.x =  element_blank(),
              panel.grid.minor.x =  element_blank(),
              panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
        xlab("Betweenness") + 
        ylab("Variables")  

# Compute clustering coefficients

clust <- clustcoef_auto(php_qnet)
clusZ.sp <- NULL
clusZ.sp$node <- rownames(clust)
clusZ.sp$clustZhang <- clust$clustZhang
clusZ.sp <- as.data.frame(clusZ.sp)

p3 <- clusZ.sp %>% 
        ggplot(aes(x= clustZhang, y = reorder(node, clustZhang))) +
        geom_point(size = 3, color ="yellow") + 
        theme_bw() +
        theme(panel.grid.major.x =  element_blank(),
              panel.grid.minor.x =  element_blank(),
              panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
        xlab("Clustering coefficient") + 
        ylab("Variables")

plot_grid(p1, p2, p3, labels=c("A", "B", "C"), ncol = 3, nrow = 1 )

abbre <- c("DH", "WH", "SS", "WM", "LF", "DEF", "BPH", "WPH", "AM", "RB", "RBB", "GLH", "STB", "BLB",  "LB", "BS", "BLS", "NBS", "RS", "LS", "SHB", "SHR", "SR", "FSM", "NB", "DP", "RTD", "RSD", "RT")
names(idn) <- abbre
InsectInjuiry <- c(1:13, 29)
Disease <- 14:28
injuries <- list(InsectInjuiry, Disease)
#==========#
cormat <- cor(idn, method = "spearman", use = "pairwise") # spearman correlation
idn_qnet <- qgraph(cormat, 
               layout = "spring", 
               threshold = 0.20,
               maximum = 1,
               groups = injuries, 
               color = c("skyblue", "wheat"),
               vsize = 5, 
               line = 3,
               diag = FALSE,
               posCol = "forestgreen",
               negCol = "firebrick3",
               borders = FALSE,
               # legend = TRUE,
               vTrans = 200,
               #nodeNames = Names,
               legend.cex = 0.3,
               title = "Indonesia"
               #filetype = "png",
               #filenames = "figs/ind_network.png"
)

qgraph(cormat, 
               layout = "spring", 
               threshold = 0.20,
               maximum = 1,
               groups = injuries, 
               color = c("skyblue", "wheat"),
               vsize = c(1.5, 10), 
               line = 3,
               diag = FALSE,
               posCol = "forestgreen",
               negCol = "firebrick3",
               borders = FALSE,
               # legend = TRUE,
               vTrans = 200,
               #nodeNames = Names,
               legend.cex = 0.3,
               title = "Indonesia"
               #filetype = "png",
               #filenames = "figs/ind_network.png"
)

# compute the topological properties
cen <- centrality_auto(idn_qnet)
topo <- cen$node.centrality
topo$node <- row.names(topo)
row.names(topo) <- NULL
topo.melt <- melt(topo)

## Using node as id variables

p1 <- topo %>% ggplot(aes(x= Strength, y = reorder(node, Strength))) +
         geom_point(size = 3, color ="red") +
        theme_bw() +
        theme(panel.grid.major.x =  element_blank(),
              panel.grid.minor.x =  element_blank(),
              panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
        xlab("Node degree") + 
        ylab("Variables")  

p2 <- topo %>% ggplot(aes(x= Betweenness, y = reorder(node, Betweenness))) + 
        geom_point(size = 3, color ="blue") +
        theme_bw() +
        theme(panel.grid.major.x =  element_blank(),
              panel.grid.minor.x =  element_blank(),
              panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
        xlab("Betweenness") + 
        ylab("Variables")  

# Compute clustering coefficients

clust <- clustcoef_auto(idn_qnet)
clusZ.sp <- NULL
clusZ.sp$node <- rownames(clust)
clusZ.sp$clustZhang <- clust$clustZhang
clusZ.sp <- as.data.frame(clusZ.sp)

p3 <- clusZ.sp %>% 
        ggplot(aes(x= clustZhang, y = reorder(node, clustZhang))) +
        geom_point(size = 3, color ="yellow") + 
        theme_bw() +
        theme(panel.grid.major.x =  element_blank(),
              panel.grid.minor.x =  element_blank(),
              panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
        xlab("Clustering coefficient") + 
        ylab("Variables")

plot_grid(p1, p2, p3, labels=c("A", "B", "C"), ncol = 3, nrow = 1 )

abbre <- c("DH", "WH", "WM", "LF", "DEF", "BPH", "WPH", "RB", "RBB", "GLH", "BLB",  "LB", "SHB", "FSM", "NB", "RT")
names(ind) <-abbre
InsectInjuiry <-c(1:10, 16)
Disease <- 11:15
injuries <- list(InsectInjuiry, Disease)

#==========#
cormat <- cor(ind, method = "spearman", use = "pairwise") # spearman correlation
ind_qnet <- qgraph(cormat, 
               layout = "spring", 
               threshold = 0.20,
               maximum = 1,
               groups = injuries, 
               color = c("skyblue", "wheat"),
               vsize = 5, 
               line = 3,
               diag = FALSE,
               posCol = "forestgreen",
               negCol = "firebrick3",
               borders = FALSE,
               # legend = TRUE,
               vTrans = 200,
               #nodeNames = Names,
               legend.cex = 0.3,
               title = "India"
               #filetype = "png",
               #filenames = "figs/ind_network.png"
)

qgraph(cormat, 
               layout = "spring", 
               threshold = 0.20,
               maximum = 1,
               groups = injuries, 
               color = c("skyblue", "wheat"),
               vsize = c(1.5, 10), 
               line = 3,
               diag = FALSE,
               posCol = "forestgreen",
               negCol = "firebrick3",
               borders = FALSE,
               # legend = TRUE,
               vTrans = 200,
               #nodeNames = Names,
               legend.cex = 0.3
               #title = "Spearman's correlation"
               #filetype = "png",
               #filenames = "figs/ind_network.png"
)

# compute the topological properties
cen <- centrality_auto(ind_qnet)
topo <- cen$node.centrality
topo$node <- row.names(topo)
row.names(topo) <- NULL
topo.melt <- melt(topo)

## Using node as id variables

p1 <- topo %>% ggplot(aes(x= Strength, y = reorder(node, Strength))) +
         geom_point(size = 3, color ="red") +
        theme_bw() +
        theme(panel.grid.major.x =  element_blank(),
              panel.grid.minor.x =  element_blank(),
              panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
        xlab("Node degree") + 
        ylab("Variables")  

p2 <- topo %>% ggplot(aes(x= Betweenness, y = reorder(node, Betweenness))) + 
        geom_point(size = 3, color ="blue") +
        theme_bw() +
        theme(panel.grid.major.x =  element_blank(),
              panel.grid.minor.x =  element_blank(),
              panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
        xlab("Betweenness") + 
        ylab("Variables")  

# Compute clustering coefficients

clust <- clustcoef_auto(ind_qnet)
clusZ.sp <- NULL
clusZ.sp$node <- rownames(clust)
clusZ.sp$clustZhang <- clust$clustZhang
clusZ.sp <- as.data.frame(clusZ.sp)

p3 <- clusZ.sp %>% 
        ggplot(aes(x= clustZhang, y = reorder(node, clustZhang))) +
        geom_point(size = 3, color ="yellow") + 
        theme_bw() +
        theme(panel.grid.major.x =  element_blank(),
              panel.grid.minor.x =  element_blank(),
              panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
        xlab("Clustering coefficient") + 
        ylab("Variables")

plot_grid(p1, p2, p3, labels=c("A", "B", "C"), ncol = 3, nrow = 1 )

abbre <- c("DH", "WH", "SS", "WM", "LF", "DEF", "BPH", "WPH", "RB", "GLH", "BLB", "LB", "BS", "BLS", "NBS", "RS", "LS", "SHB", "SHR", "SR", "FSM", "NB", "DP", "RT")
names(tha) <- abbre
InsectInjuiry <- c(1:10, 24)
Disease <- 11:23
injuries <- list(InsectInjuiry, Disease)

#==========#
cormat <- cor(tha, method = "spearman", use = "pairwise") # spearman correlation
tha_qnet <- qgraph(cormat, 
               layout = "spring", 
               threshold = 0.20,
               maximum = 1,
               groups = injuries, 
               color = c("skyblue", "wheat"),
               vsize = 5, 
               line = 3,
               diag = FALSE,
               posCol = "forestgreen",
               negCol = "firebrick3",
               borders = FALSE,
               # legend = TRUE,
               vTrans = 200,
               #nodeNames = Names,
               legend.cex = 0.3,
               title = "Thailand"
               #filetype = "png",
               #filenames = "figs/ind_network.png"
)

qgraph(cormat, 
               layout = "spring", 
               threshold = 0.20,
               maximum = 1,
               groups = injuries, 
               color = c("skyblue", "wheat"),
               vsize = c(1.5, 10), 
               line = 3,
               diag = FALSE,
               posCol = "forestgreen",
               negCol = "firebrick3",
               borders = FALSE,
               # legend = TRUE,
               vTrans = 200,
               #nodeNames = Names,
               legend.cex = 0.3,
               title = "Thailand"
               #filetype = "png",
               #filenames = "figs/ind_network.png"
)

# compute the topological properties
cen <- centrality_auto(tha_qnet)
topo <- cen$node.centrality
topo$node <- row.names(topo)
row.names(topo) <- NULL
topo.melt <- melt(topo)

## Using node as id variables

p1 <- topo %>% ggplot(aes(x= Strength, y = reorder(node, Strength))) +
         geom_point(size = 3, color ="red") +
        theme_bw() +
        theme(panel.grid.major.x =  element_blank(),
              panel.grid.minor.x =  element_blank(),
              panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
        xlab("Node degree") + 
        ylab("Variables")  

p2 <- topo %>% ggplot(aes(x= Betweenness, y = reorder(node, Betweenness))) + 
        geom_point(size = 3, color ="blue") +
        theme_bw() +
        theme(panel.grid.major.x =  element_blank(),
              panel.grid.minor.x =  element_blank(),
              panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
        xlab("Betweenness") + 
        ylab("Variables")  

# Compute clustering coefficients

clust <- clustcoef_auto(tha_qnet)
clusZ.sp <- NULL
clusZ.sp$node <- rownames(clust)
clusZ.sp$clustZhang <- clust$clustZhang
clusZ.sp <- as.data.frame(clusZ.sp)

p3 <- clusZ.sp %>% 
        ggplot(aes(x= clustZhang, y = reorder(node, clustZhang))) +
        geom_point(size = 3, color ="yellow") + 
        theme_bw() +
        theme(panel.grid.major.x =  element_blank(),
              panel.grid.minor.x =  element_blank(),
              panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
        xlab("Clustering coefficient") + 
        ylab("Variables")

plot_grid(p1, p2, p3, labels=c("A", "B", "C"), ncol = 3, nrow = 1 )

abbre <- c("DH", "WH", "SS", "WM", "LF", "DEF", "BPH", "WPH", "AM", "RB", "GLH", "BLB", "LB", "BS", "BLS", "NBS", "RS", "LS", "SHB", "SHR", "SR", "FSM", "NB", "DP", "RTD", "RDS", "RT")
names(vnm) <-abbre
InsectInjuiry <- c(1:11, 27)
Disease <- 12:26
injuries <- list(InsectInjuiry, Disease)

#==========#
cormat <- cor(vnm, method = "spearman", use = "pairwise") # spearman correlation
vnm_qnet <- qgraph(cormat, 
               layout = "spring", 
               threshold = 0.20,
               maximum = 1,
               groups = injuries, 
               color = c("skyblue", "wheat"),
               vsize = 5, 
               line = 3,
               diag = FALSE,
               posCol = "forestgreen",
               negCol = "firebrick3",
               borders = FALSE,
               # legend = TRUE,
               vTrans = 200,
               #nodeNames = Names,
               legend.cex = 0.3,
               title = "Vietnam"
               #filetype = "png",
               #filenames = "figs/ind_network.png"
)

qgraph(cormat, 
               layout = "spring", 
               threshold = 0.20,
               maximum = 1,
               groups = injuries, 
               color = c("skyblue", "wheat"),
               vsize = c(1.5, 10), 
               line = 3,
               diag = FALSE,
               posCol = "forestgreen",
               negCol = "firebrick3",
               borders = FALSE,
               # legend = TRUE,
               vTrans = 200,
               #nodeNames = Names,
               legend.cex = 0.3,
               title = "Vietnam"
               #filetype = "png",
               #filenames = "figs/ind_network.png"
)

# compute the topological properties
cen <- centrality_auto(vnm_qnet)
topo <- cen$node.centrality
topo$node <- row.names(topo)
row.names(topo) <- NULL
topo.melt <- melt(topo)

## Using node as id variables

p1 <- topo %>% ggplot(aes(x= Strength, y = reorder(node, Strength))) +
         geom_point(size = 3, color ="red") +
        theme_bw() +
        theme(panel.grid.major.x =  element_blank(),
              panel.grid.minor.x =  element_blank(),
              panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
        xlab("Node degree") + 
        ylab("Variables")  

p2 <- topo %>% ggplot(aes(x= Betweenness, y = reorder(node, Betweenness))) + 
        geom_point(size = 3, color ="blue") +
        theme_bw() +
        theme(panel.grid.major.x =  element_blank(),
              panel.grid.minor.x =  element_blank(),
              panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
        xlab("Betweenness") + 
        ylab("Variables")  

# Compute clustering coefficients

clust <- clustcoef_auto(vnm_qnet)
clusZ.sp <- NULL
clusZ.sp$node <- rownames(clust)
clusZ.sp$clustZhang <- clust$clustZhang
clusZ.sp <- as.data.frame(clusZ.sp)

p3 <- clusZ.sp %>% 
        ggplot(aes(x= clustZhang, y = reorder(node, clustZhang))) +
        geom_point(size = 3, color ="yellow") + 
        theme_bw() +
        theme(panel.grid.major.x =  element_blank(),
              panel.grid.minor.x =  element_blank(),
              panel.grid.major.y = element_line(color = "grey", linetype = 3)) +
        xlab("Clustering coefficient") + 
        ylab("Variables")

plot_grid(p1, p2, p3, labels=c("A", "B", "C"), ncol = 3, nrow = 1 )